"""
    DiagnosticFileFormat

GSI diagnostic file format specifications and data structures.

This module documents and implements the binary format used by NOAA/NCEP GSI
for diagnostic output files. The format is based on Fortran sequential
unformatted binary files.

# File Types

## Observation Space Diagnostics
- `diag_obs_setup_anl.YYYYMMDDHH` - Analysis diagnostics
- `diag_obs_setup_ges.ensmean` - Background/guess diagnostics
- `diag_[obstype]_anl.[iteration]` - Per-type analysis diagnostics
- `diag_[obstype]_ges.[iteration]` - Per-type background diagnostics

## Statistics Files (fort.*)
- `fort.201` - Surface pressure fit statistics
- `fort.202` - Wind fit statistics
- `fort.203` - Temperature fit statistics
- `fort.204` - Moisture fit statistics
- `fort.205-237` - Additional observation type statistics

# Binary Format Structure

## Header Record
```
Record 1: Header (38 bytes)
  - Variable name (10 bytes, char)
  - Variable type (10 bytes, char)
  - Padding/metadata (18 bytes)
```

## Grid Information Record
```
Record 2: Grid metadata (variable size)
  - Grid dimensions (nx, ny, nz)
  - Grid type code
  - Coordinate arrays
```

## Observation Records
```
Records 3+: Individual observation data
Each record contains multiple observations of the same type
  - Latitude (Float64)
  - Longitude (Float64)
  - Pressure/height (Float64)
  - Observation value (Float64)
  - Observation error (Float64)
  - Background value H(xb) (Float64)
  - Analysis value H(xa) (Float64)
  - Quality control flag (Int32)
  - Use flag (Int32)
  - Observation type code (Int32)
  - Station ID or channel (16 bytes, char)
  - Time offset (Float64)
  - Additional metadata...
```
"""
module DiagnosticFileFormat

using Printf
using Dates

export DiagnosticHeader, ObservationRecord, GridMetadata
export create_diagnostic_header, pack_observation_record, unpack_observation_record
export OBSERVATION_TYPES, OBSERVATION_SUBTYPES

# GSI Observation Type Codes
const OBSERVATION_TYPES = Dict{String, Int32}(
    "ps" => 120,        # Surface pressure
    "t" => 130,         # Temperature
    "q" => 133,         # Specific humidity
    "u" => 131,         # U-wind
    "v" => 132,         # V-wind
    "gps" => 153,       # GPS refractivity
    "rad" => 175,       # Radiance
    "conv" => 180,      # Conventional data
    "light" => 195,     # Lightning
)

const OBSERVATION_SUBTYPES = Dict{String, Int32}(
    "adpupa" => 120,    # Radiosonde
    "aircar" => 130,    # Aircraft
    "aircft" => 131,    # Aircraft
    "satwnd" => 240,    # Satellite winds
    "profl" => 220,     # Wind profiler
    "vadwnd" => 224,    # VAD winds
    "satemp" => 180,    # Satellite retrieval
    "adpsfc" => 181,    # Surface obs
    "sfcshp" => 180,    # Ship obs
    "sfcbog" => 187,    # Bogus surface
    "spssmi" => 189,    # SSM/I
    "syndat" => 190,    # Synthetic
    "goesnd" => 156,    # GOES sounder
    "gpsref" => 153,    # GPS refractivity
)

"""
    DiagnosticHeader

Header record for GSI diagnostic file.

# Fields
- `var_name`: Variable name (e.g., "ps", "t", "q")
- `var_type`: Variable type descriptor
- `nx`: Grid dimension x
- `ny`: Grid dimension y
- `nz`: Grid dimension z
- `analysis_time`: Analysis date/time
- `n_observations`: Number of observations
"""
struct DiagnosticHeader
    var_name::String
    var_type::String
    nx::Int32
    ny::Int32
    nz::Int32
    analysis_time::DateTime
    n_observations::Int
end

"""
    GridMetadata

Grid metadata for diagnostic file.

# Fields
- `grid_type`: Grid type code (1=regional, 2=global lat-lon, 3=spectral)
- `nx`: Number of longitude points
- `ny`: Number of latitude points
- `nz`: Number of vertical levels
- `dx`: Grid spacing in x (degrees or meters)
- `dy`: Grid spacing in y (degrees or meters)
- `lon_min`: Minimum longitude
- `lat_min`: Minimum latitude
- `vertical_coords`: Vertical coordinate values (pressure or sigma)
"""
struct GridMetadata
    grid_type::Int32
    nx::Int32
    ny::Int32
    nz::Int32
    dx::Float64
    dy::Float64
    lon_min::Float64
    lat_min::Float64
    vertical_coords::Vector{Float64}
end

"""
    ObservationRecord

Single observation record in GSI diagnostic file.

# Fields
- `lat`: Latitude (degrees)
- `lon`: Longitude (degrees)
- `pressure`: Pressure level (hPa) or height (m)
- `obs_value`: Observed value
- `obs_error`: Observation error standard deviation
- `background_value`: Background (first guess) value at obs location H(xb)
- `analysis_value`: Analysis value at obs location H(xa)
- `qc_flag`: Quality control flag (0=good, >0=rejected)
- `use_flag`: Use in analysis flag (1=used, 0=not used, -1=monitored)
- `obs_type`: Observation type code
- `obs_subtype`: Observation subtype code
- `station_id`: Station ID or satellite channel
- `time_offset`: Time offset from analysis time (hours)
- `ob_minus_background`: Innovation (O-B)
- `ob_minus_analysis`: Analysis residual (O-A)
- `inverse_obs_error`: Inverse observation error (1/σ)
- `variational_qc_weight`: Variational QC weight (0-1)
"""
mutable struct ObservationRecord
    lat::Float64
    lon::Float64
    pressure::Float64
    obs_value::Float64
    obs_error::Float64
    background_value::Float64
    analysis_value::Float64
    qc_flag::Int32
    use_flag::Int32
    obs_type::Int32
    obs_subtype::Int32
    station_id::String
    time_offset::Float64
    ob_minus_background::Float64
    ob_minus_analysis::Float64
    inverse_obs_error::Float64
    variational_qc_weight::Float64
end

"""
    ObservationRecord()

Create an empty observation record with default values.
"""
function ObservationRecord()
    return ObservationRecord(
        0.0,    # lat
        0.0,    # lon
        0.0,    # pressure
        0.0,    # obs_value
        1.0,    # obs_error
        0.0,    # background_value
        0.0,    # analysis_value
        0,      # qc_flag
        0,      # use_flag
        0,      # obs_type
        0,      # obs_subtype
        "",     # station_id
        0.0,    # time_offset
        0.0,    # ob_minus_background
        0.0,    # ob_minus_analysis
        1.0,    # inverse_obs_error
        1.0     # variational_qc_weight
    )
end

"""
    create_diagnostic_header(var_name::String, grid_size::Tuple,
                            analysis_time::DateTime, n_obs::Int)

Create a diagnostic file header.
"""
function create_diagnostic_header(var_name::String, grid_size::Tuple,
                                  analysis_time::DateTime, n_obs::Int)
    nx, ny, nz = length(grid_size) == 3 ? grid_size : (grid_size..., 1)

    return DiagnosticHeader(
        var_name,
        var_name,  # var_type same as var_name for simplicity
        Int32(nx),
        Int32(ny),
        Int32(nz),
        analysis_time,
        n_obs
    )
end

"""
    pack_observation_record(obs::ObservationRecord)

Pack observation record into binary format (Vector{UInt8}).

Returns a byte array suitable for writing to Fortran binary file.
The format matches GSI's observation diagnostic record structure.
"""
function pack_observation_record(obs::ObservationRecord)
    # Calculate total record size
    # Standard GSI observation record:
    # - 17 Float64 values (8 bytes each) = 136 bytes
    # - 3 Int32 values (4 bytes each) = 12 bytes
    # - 1 String (16 bytes) = 16 bytes
    # Total = 164 bytes (approximate, may vary by observation type)

    io = IOBuffer()

    # Write floating point fields (Float64, big-endian)
    write(io, hton(obs.lat))
    write(io, hton(obs.lon))
    write(io, hton(obs.pressure))
    write(io, hton(obs.obs_value))
    write(io, hton(obs.obs_error))
    write(io, hton(obs.background_value))
    write(io, hton(obs.analysis_value))
    write(io, hton(obs.time_offset))
    write(io, hton(obs.ob_minus_background))
    write(io, hton(obs.ob_minus_analysis))
    write(io, hton(obs.inverse_obs_error))
    write(io, hton(obs.variational_qc_weight))

    # Write integer fields (Int32, big-endian)
    write(io, hton(obs.qc_flag))
    write(io, hton(obs.use_flag))
    write(io, hton(obs.obs_type))
    write(io, hton(obs.obs_subtype))

    # Write station ID (fixed 16-byte string, space-padded)
    station_id_padded = rpad(obs.station_id, 16)[1:16]
    write(io, Vector{UInt8}(station_id_padded))

    return take!(io)
end

"""
    unpack_observation_record(data::Vector{UInt8})

Unpack observation record from binary format.

Returns an ObservationRecord struct populated from the byte array.
"""
function unpack_observation_record(data::Vector{UInt8})
    io = IOBuffer(data)

    obs = ObservationRecord()

    # Read floating point fields (Float64, big-endian)
    obs.lat = ntoh(read(io, Float64))
    obs.lon = ntoh(read(io, Float64))
    obs.pressure = ntoh(read(io, Float64))
    obs.obs_value = ntoh(read(io, Float64))
    obs.obs_error = ntoh(read(io, Float64))
    obs.background_value = ntoh(read(io, Float64))
    obs.analysis_value = ntoh(read(io, Float64))
    obs.time_offset = ntoh(read(io, Float64))
    obs.ob_minus_background = ntoh(read(io, Float64))
    obs.ob_minus_analysis = ntoh(read(io, Float64))
    obs.inverse_obs_error = ntoh(read(io, Float64))
    obs.variational_qc_weight = ntoh(read(io, Float64))

    # Read integer fields (Int32, big-endian)
    obs.qc_flag = ntoh(read(io, Int32))
    obs.use_flag = ntoh(read(io, Int32))
    obs.obs_type = ntoh(read(io, Int32))
    obs.obs_subtype = ntoh(read(io, Int32))

    # Read station ID (fixed 16-byte string)
    station_bytes = read(io, 16)
    obs.station_id = String(station_bytes) |> rstrip

    return obs
end

"""
    get_obs_type_code(obs_type_name::String)

Get observation type code from name.
"""
function get_obs_type_code(obs_type_name::String)
    return get(OBSERVATION_TYPES, lowercase(obs_type_name), Int32(0))
end

"""
    get_obs_subtype_code(obs_subtype_name::String)

Get observation subtype code from name.
"""
function get_obs_subtype_code(obs_subtype_name::String)
    return get(OBSERVATION_SUBTYPES, lowercase(obs_subtype_name), Int32(0))
end

"""
    compute_innovation_statistics(observations::Vector{ObservationRecord})

Compute innovation statistics from observation records.

Returns a dictionary with:
- "mean_ob": Mean of observations
- "mean_bg": Mean of background values
- "mean_innovation": Mean innovation (O-B)
- "rms_innovation": RMS innovation
- "bias": Bias (same as mean innovation)
- "std_innovation": Standard deviation of innovations
- "n_used": Number of observations used
- "n_rejected": Number of observations rejected
"""
function compute_innovation_statistics(observations::Vector{ObservationRecord})
    n_total = length(observations)
    n_used = count(obs -> obs.use_flag > 0, observations)
    n_rejected = count(obs -> obs.qc_flag > 0, observations)

    # Get used observations only
    used_obs = filter(obs -> obs.use_flag > 0, observations)

    if isempty(used_obs)
        return Dict(
            "mean_ob" => 0.0,
            "mean_bg" => 0.0,
            "mean_innovation" => 0.0,
            "rms_innovation" => 0.0,
            "bias" => 0.0,
            "std_innovation" => 0.0,
            "n_total" => n_total,
            "n_used" => 0,
            "n_rejected" => n_rejected
        )
    end

    innovations = [obs.ob_minus_background for obs in used_obs]
    obs_values = [obs.obs_value for obs in used_obs]
    bg_values = [obs.background_value for obs in used_obs]

    mean_innovation = sum(innovations) / length(innovations)
    rms_innovation = sqrt(sum(innovations.^2) / length(innovations))
    std_innovation = sqrt(sum((innovations .- mean_innovation).^2) / length(innovations))

    return Dict(
        "mean_ob" => sum(obs_values) / length(obs_values),
        "mean_bg" => sum(bg_values) / length(bg_values),
        "mean_innovation" => mean_innovation,
        "rms_innovation" => rms_innovation,
        "bias" => mean_innovation,
        "std_innovation" => std_innovation,
        "n_total" => n_total,
        "n_used" => n_used,
        "n_rejected" => n_rejected
    )
end

end # module DiagnosticFileFormat
